**Better grid search for \alpha
**Keep only preferred scobit

qui {

clear all

program define MESim, rclass

	syntax [, b0(real 1) b1(real 1) b2(real 1) b3(real 1) theta0(real 1) theta1(real 1) alpha(real 1) nobs(integer 10) t(integer 10)]
	
	clear
	set obs `nobs' 
	ret sca nobs = `nobs'
	ret sca t = `t'
	ret sca beta0 = `b0'
	ret sca theta0 = `theta0'
	ret sca theta1 = `theta1'
	ret sca alpha = `alpha'
	
	g id = _n
	g fe = rnormal(0,5)
	expand `t'
	bys id: g t=_n
	xtset i t
	
	g z 	= rpoisson(9)
	g x1	= rchi2(25)/100				// slope
	g x2	= rchi2(45)					// dist_any_rooad
	
	g u 	= rnormal(0,5)
	g d 	= (-8+.1*z-.1*x1+0.05*x2+0.5*fe+u>0)
	noi tab d
	noi corr d x1 x2 z
	
	g mu 	= runiform()
	g p 	= 1 - (1/(1+exp(`b0' + `b1'*x1 + `b2'*x2 + `b3'*d + fe))^`alpha')
	g ystar = (mu<p) 	// truth, binary  
	
	g m 	= `b0' + `b1'*x1 + `b2'*x2 + `b3'*d + fe
	g m1 	= `b0' + `b1'*x1 + `b2'*x2 + `b3' + fe
	g m0 	= `b0' + `b1'*x1 + `b2'*x2 + fe
	g me_d0 = ( 1 - (1/(1+exp(m1))^`alpha')) - ( 1 - (1/(1+exp(m0))^`alpha'))  
	g me_x10 = `b1'*exp(m)*`alpha'*((1+exp(m))^((-1)*`alpha' - 1))
	g me_x20 = `b2'*exp(m)*`alpha'*((1+exp(m))^((-1)*`alpha' - 1))
	drop m1 m0	
	
	g u1	= uniform()
	g u2	= uniform()
	g y 	= ystar
	g fp 	= normal(`theta0'-0.10*z)
	g fn 	= normal(`theta1'+0.15*z)
	replace y = 1 if ystar==0 & u1<fp			// mismeasured, binary
	replace y = 0 if ystar==1 & u2<fn			// mismeasured, binary

	tab ystar
	count if ystar==0
	loc n0 = r(N)
	count if ystar==1
	loc n1 = r(N)
	count if y==1 & ystar==0
	ret sca fp = r(N)/`n0'
	count if y==0 & ystar==1
	ret sca fn = r(N)/`n1'
		
	**** CRE CONTROLS
	bys id: egen mx1 = mean(x1)
	bys id: egen mx2 = mean(x2)
	bys id: egen md = mean(d)
	bys id: egen mz = mean(z)		
		
	**** BINARY TRUE (CRE)
	logit ystar i.d x1 x2 md mx1 mx2, nolog 
	cap drop xb
	predict xb, xb
	g me_d1 = logistic(xb+_b[1.d]*(1-d)) - logistic(xb-_b[1.d]*d) 				
	g me_x11 = logisticden(xb)*_b[x1]
	g me_x21 = logisticden(xb)*_b[x2]
	
	**** BINARY MISMEASURED
	logit y i.d x1 x2 md mx1 mx2, nolog 
	loc a = _b[1.d]
	loc b1 = _b[x1]
	loc b2 = _b[x2]
	loc c = _b[_cons]
	loc m0 = _b[md]
	loc m1 = _b[mx1]
	loc m2 = _b[mx2]
	loc from "y:_cons=`c' y:1.d=`a' y:x1=`b1' y:x2=`b2' y:md=`m0' y:mx1=`m1' y:mx2=`m2'"	
	cap drop xb
	predict xb, xb
	g me_d2 = logistic(xb+_b[1.d]*(1-d)) - logistic(xb-_b[1.d]*d) 				
	g me_x12 = logisticden(xb)*_b[x1]
	g me_x22 = logisticden(xb)*_b[x2]
	
	**** BINARY AD HOC SOLUTION
	logit y i.d x1 x2 z md mx1 mx2 mz, nolog 
	cap drop xb
	predict xb, xb
	g me_d3 = logistic(xb+_b[1.d]*(1-d)) - logistic(xb-_b[1.d]*d) 				
	g me_x13 = logisticden(xb)*_b[x1]
	g me_x23 = logisticden(xb)*_b[x2]
		
	**** BINARY SOLUTION
	cap mclogit y i.d x1 x2 md mx1 mx2, mc(z) diff from(a0:_cons=-1 a0:z=-.1 a1:_cons=-1 a1:z=.1 `from') tech(nr 5 dfp 10 bhhh 10) 
	if e(converged) == 1 {
		cap drop xb
		predict xb, xb
		g me_d4 = logistic(xb+_b[1.d]*(1-d)) - logistic(xb-_b[1.d]*d) 				
		g me_x14 = logisticden(xb)*_b[x1]
		g me_x24 = logisticden(xb)*_b[x2]
	}
	else {
		g me_d4 = .
		g me_x14 = .
		g me_x24 = .
	}
				
	**** LPM MISMEASURED
	xtreg y d x1 x2, fe
	g me_d5 = _b[d]
	g me_x15 = _b[x1]
	g me_x25 = _b[x2]
	
	**** LPM AD HOC SOLUTION
	xtreg y d x1 x2 z, fe
	g me_d6 = _b[d]
	g me_x16 = _b[x1]
	g me_x26 = _b[x2]
	

	
	
	
	
	
	**** MC SCOBIT
	loc ea = ln(.5)
	constraint 3 /lnalpha = `ea'
	constraint 4 [ln_a2]_cons = `ea'
	loc ea = ln(.75)
	constraint 5 /lnalpha = `ea'
	constraint 6 [ln_a2]_cons = `ea'
	
	**** MC SCOBIT 
	loc logL = -999999999
	g me_d7 = .
	g me_x17 = .
	g me_x27 = .
	forval gs = 0.2(0.05)0.95 {
	
		loc ea = ln(`gs')
		constraint 1 /lnalpha = `ea'
		constraint 2 [ln_a2]_cons = `ea'
		scobit y i.d x1 x2 md mx1 mx2, nolog constraint(1) 
		loc a = _b[1.d]
		loc b1 = _b[x1]
		loc b2 = _b[x2]
		loc c = _b[_cons]
		loc m0 = _b[md]
		loc m1 = _b[mx1]
		loc m2 = _b[mx2]
		loc from "y:_cons=`c' y:1.d=`a' y:x1=`b1' y:x2=`b2' y:md=`m0' y:mx1=`m1' y:mx2=`m2'"	
							
		cap mcre y i.d x1 x2 md mx1 mx2, mc(z) diff from(a0:_cons=-1 a0:z=-.1 a1:_cons=-1 a1:z=.1 `from') constraint(2) tech(nr 5 dfp 10 bhhh 10)	
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc alp = exp([ln_a2]_cons)	
				replace me_d7 = 1 - (1/(1+exp(xb+_b[1.d]*(1-d)))^`alp') - ( 1 - (1/(1+exp(xb-_b[1.d]*d))^`alp') )
				replace me_x17 = [y]x1*exp(xb)*`alp'*((1+exp(xb))^((-1)*`alp' - 1))
				replace me_x27 = [y]x2*exp(xb)*`alp'*((1+exp(xb))^((-1)*`alp' - 1))
				loc logL = e(ll)
			}
		}
		
	}
	
	**** RETURN RESULTS
	forval j = 0(1)7 {
		su me_d`j', meanonly
		ret sca ame_d_`j' = r(mean)
		su me_x1`j', meanonly
		ret sca ame_x1_`j' = r(mean)
		su me_x2`j', meanonly
		ret sca ame_x2_`j' = r(mean)
	}
	
end //end program define

#delimit ;
	
simulate
	nobs   		= r(nobs)
	T			= r(t)
	alpha		= r(alpha)
	beta0		= r(beta0)
	theta0		= r(theta0)
	theta1		= r(theta1)
	fp			= r(fp)
	fn			= r(fn)
	ame_d_0		= r(ame_d_0)
	ame_x1_0	= r(ame_x1_0)
	ame_x2_0	= r(ame_x2_0)
	ame_d_1		= r(ame_d_1)
	ame_x1_1	= r(ame_x1_1)
	ame_x2_1	= r(ame_x2_1)
	ame_d_2		= r(ame_d_2)
	ame_x1_2	= r(ame_x1_2)
	ame_x2_2	= r(ame_x2_2)
	ame_d_3		= r(ame_d_3)
	ame_x1_3	= r(ame_x1_3)
	ame_x2_3	= r(ame_x2_3)
	ame_d_4		= r(ame_d_4)
	ame_x1_4	= r(ame_x1_4)
	ame_x2_4	= r(ame_x2_4)
	ame_d_5		= r(ame_d_5)
	ame_x1_5	= r(ame_x1_5)
	ame_x2_5	= r(ame_x2_5)
	ame_d_6		= r(ame_d_6)
	ame_x1_6	= r(ame_x1_6)
	ame_x2_6	= r(ame_x2_6)
	ame_d_7		= r(ame_d_7)
	ame_x1_7	= r(ame_x1_7)
	ame_x2_7	= r(ame_x2_7)
	,		
	seed(${seed}) reps(${reps}) saving(MEBinaryDGP${dgp}.dta, replace) nodots : MESim, 
	b0(${beta0}) b1(${beta1}) b2(${beta2}) b3(${beta3}) theta0(${theta0}) theta1(${theta1}) 
	alpha(${alpha}) nobs(${nobs}) t(${T});
			
#delimit cr


} 
//end qui
